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^ ■ Abstract 
O 

o, 

\ We have simulated lattice QCD directly in the chiral limit of zero quark mass 

I by adding an additional, irrelevant 4-fermion interaction to the standard ac- 

& , tion. Using lattices having temporal extent of six and spatial extents of twelve 

Qh! and eighteen, we find that the theory with 2 massless staggered quark flavors 



has a second order finite temperature phase transition. The critical expo- 



nents Pmag, 8 an d v are measured and favour tricritical behaviour over that 

h : 

expected by universality arguments. The pion screening mass is consistent 

with zero below the transition, but is degenerate with the nonzero <r(/o) mass 

above the transition, indicating the restoration of chiral symmetry. 
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There are several bottlenecks in the simulation studies of lattice field theories by 
standard algorithms [I], 0. One of the most harmful is that the chiral limit cannot be 
simulated directly because the standard algorithms fail to converge in the limit of vanishing 
bare fermion mass and become prohibitively expensive at small quark masses. Various chiral 
transitions, such as the finite temperature Quantum Chromodynamics (QCD) transition to 
be discussed here, have been difficult to study quantitatively because of this. In particu- 
lar, studies of critical scaling and the extraction of critical exponents have been especially 
difficult. This has led people to extract the critical exponents associated with the mass de- 
pendence of the heights and positions of susceptibility peaks [§,|J and to compare the scaling 
behaviour of the chiral condensate as both the coupling constant and the quark mass are 
varied ||, rather than studying the scaling of order parameters in terms of the coupling 
constant and mass independently. The most extensive data has been obtained for N t = 4 
where it differs significantly from what universality arguments predict. 

These problems are largely solved by adding a four-fermi interaction to standard 
lattice actions with staggered fermions ||. The resulting lattice action can be simulated 
directly in the chiral limit because an auxiliary scalar field a (essentially the chiral condensate 
(^i/j)) acts as a dynamical ('constituent') quark mass term insuring that the inversion of 
the Dirac operator will be successful and very fast. This approach is far more efficient and 
physical than the traditional lattice algorithm in which the chiral limit is singular and a 
nonzero bare quark mass m is a necessity. 

Using this strategy we have simulated lattice Quantum Chromodynamics (QCD) 
with two flavors of staggered quarks on 12 3 x 6 and 18 3 x 6 lattices, at zero quark mass, in 
order to determine the position and nature of the finite temperature transition. In addition 
to measuring the standard order parameters, we have measured the pion and cr(/ ) screening 
masses to probe the nature of chiral symmetry restoration at this transition. We find that 
the two flavor theory experiences a second order chiral symmetry restoring phase transition. 
We make independent determinations of the critical indices (3 magi u, and 5. (3 ma g differs 
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significantly from that of the 0(4) (or 0(2)) sigma model in 3 dimensions, expected in the 
standard scenario of dimensional reduction and universality for this transition J7|, but is in 
excellent agreement with that of the 3-dimensonal tricritical point. Since we have found that 
the N t = 4 transition is first order at 7 = 10 || and a tricritical point would be expected to 
separate a line of first order transitions from a line of second order transitions, this suggests 
that this theory would have a normal critical point with the critical exponents expected 
from universality for N t > 8. If so, the anomalous behaviour of the N t = 4,6 transitions is 
due to irrelevant terms in the lattice action which are important only on coarse lattices. 

In this letter we review the 'standard' expectations and alternative scenarios for 
the finite temperature transition in continuum and lattice QCD. We review the lattice 
formulation of QCD with chiral four-fermion interactions, and present our data and fits. 
We end with discussions and conclusions, and directions for future research. 

Dimensional reduction and universality predict that 2-flavour QCD has the critical 
behaviour of the 3-dimensional 0(4) sigma model |7]]. This theory has a second order phase 
transition with critical indices |J, v = .73(2), f3 mag = .38(1), 5 = 4.82(5), etc. The lack 
of full flavour symmetry of the staggered fermion method is expected to reduce this to an 
0(2) sigma model with critical indices v = .71(2), j3 mag = .35(1) and 5 = 4.81(1) ||. A 
breakdown of such arguments or the presence of lattice artifacts can produce a first order 
transition. If both a first and simple second order transitions are possible, so is a tricritical 
point for which v = |, (3 mag = | and 6 = 5 or 3, the second value being for a second 
independent symmetry breaking operator this more complex transition admits |10] . Finally, 
if dimensional reduction fails to occur, one would expect a mean field transition with v = |, 
Pmag = \ and S = 3. 

This suggests that a direct measurement of f3 mag , requiring simulations at m = will 
be crucial in determining the nature of the finite temperature chiral transition in 2-flavour 
(lattice) QCD. 

The lattice QCD action used in this study contains an additional chirally invariant 
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4-fermion interaction. Because the extra interaction is irrelevant, such an action should lie 
in the same universality class as the standard action. Ideally, such an interaction should be 
chosen to have the SU(Nf) x SU(Nf) flavour symmetry of the original QCD action. However, 
when one introduces auxiliary scalar and pseudoscalar fields to render this action quadratic 
in the fermion fields — which is necessary for lattice simulations — , the fermion determinant 
is no longer real, even in the continuum limit. Thus for 2-flavour QCD (Nf = 2), we make a 
simpler choice and choose a 4-fermion term with the symmetry U (1) x £7(1) C 577(2) x SU(2), 
where £7(1) x £7(1) is generated by ^3,75X3), which preserves the symmetries of the lattice 
action. The euclidean Lagrangian density for this theory is 

C = \f^ u + $(fl+ a + z7r 75 T 3 + m)ij + ^(a 2 + tt 2 ) (1) 

when we have introduced auxiliary fields a and ir to render it quadratic in the fermion fields. 
7 = 3/A 2 , where the quartic term in the original Lagrangian had coefficient X 2 /6Nf. The 
molecular dynamics Lagrangian for a particular staggered fermion lattice transcription of 
this theory has been given in reference ||. There we reported simulation results on lattices 
with N t — 4 and discussed the choice of parameters as well as systematic errors. Here we 
shall discuss higher precision simulations on larger lattices, 18 3 x 6 and 12 3 x 6. 

By doing simulations on two lattice sizes we were able to isolate finite volume effects 
in several observables and check that they had the behaviour expected of finite size scaling 
theory. We ran simulations at 7 = 10 and 20 on a 12 3 x 6 and 7 = 20 on an 18 3 x 6 lattice. 
We carefully checked that our results do not depend on 7 or on the lattice size. These 
technically important results will be discussed at greater length elsewhere ||. 

Since our m = lattice action does not single out a preferred chiral direction, we 
accumulate the magnitude of the order parameter, y (V^) 2 + (^75 r 3'0) 2 > rather than the 
order parameter itself, which averages to zero. This is not a true order parameter since it 
does not vanish identically in the symmetric phase. However, for a large system of volume V, 
its values in the high temperature phase are of order 1 / \JV . In the low temperature phase 
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the order parameter can be fit to a standard critical form, A{j3 c — f3Y maa , if (3 is chosen 
within the model's scaling window near f3 c . Of course, we cannot simulate the model with 
(3 chosen too close to /3 C where the order parameter is very small, and finite size rounding in 
this effective order parameter is considerable. This forces us to work slightly away from the 
critical point and would produce indecisive results if the scaling window were particularly 
small. These potential problems bear close scrutiny. 

The simulations of the lattice version of Equation |l| were performed using the hybrid 
molecular-dynamics algorithm with "noisy" fermions allowing us to tune Nf to 2 flavours 
||. The simulations on an 18 3 x 6 lattice were performed at 14 (3 values from 5.39 to 5.45, 
the observed chiral transition occurring close to (3 — 5.423. This enabled us to fit order 
parameters in the scaling window on both sides of the transition. Runs at single /3's varied 
in length from 5,000 to 50,000 molecular-dynamics time units close to the critical point (this 



computer "time" is the same as that used in related works by other groups, such as |TI|). We 
were forced to make such lengthy runs because of critical slowing down. Order parameters 
measured included the local chiral condensates (ipip) and (ip^Tsip) , their counterparts in 
terms of the auxiliary fields viz (a) and (it), the thermal Wilson line, the plaquette, and 
operators contributing to the partial entropies of each field. These were measured every 2 
time units and binned to minimize the effect of correlations on error estimates. The molecular 
dynamics "time" increment for updating was dt = 0.05 ||. This was large enough that the 
order parameters (a) and (n) have sizable 0(dt 2 ) errors. However, since these errors can be 
absorbed in a redefinition of the couplings in the action, they should not affect universality. 



The data for the magnitude of the order parameter y (ipip) 2 + (V'75' r 3V ; ) 2 is plotted 
against f3 in Figure |l|. The points between (3 values of 5.41 and 5.4225 were fitted to a 
simple powerlaw A(f3 c — (3) l3ma3 . The fit has a confidence level of 94 percent and determined 
the parameters (3 C = 5.4230(2), a = 1.19(12) and, most importantly, the critical index 
Pmag = 0.27(2), in good agreement with the 0.25 of the tricritical point, but inconsistent 
with the 0(4), 0(2), and mean field values. 
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FIG. 1. The order parameter plotted against [3. The solid line is the fit discussed in the text. 

We see in the figure that points at stronger coupling than (3 = 5.41 are outside the 
scaling window. For example, if the points at (3 = 5.4 and (3 = 5.405 were included in the 
fit, the confidence level would fall to 28 percent while the predicted best value for (3 mag does 
not change significantly. Our simulations on 12 3 x 6 lattices produced values for the order 
parameter which were essentially identical to those shown in the figure for (3 < 5.4225, but 
were considerably larger for (3 > 5.4225. Using the 2 lattice sizes to remove the finite volume 
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effects gave f3 mag = 0.24(2) consistent with the above. The data on the weak coupling side of 
the transition is, in fact, compatible with the finite size scaling prediction that the magnitude 
of the order parameter should fall to zero as 1/y/Vs, where V s is the volume of the spatial 
box. For the 12 3 x 6 lattice and 7 = 10, our fits predicted f3 mag = 0.27(3), in good agreement 
with the 18 3 x 6, 7 = 20 results, while f3 c = 5.4650(1). 

A second independent critical index is 5 which controls the size of the order parameter 
at the critical temperature as the quark mass is turned on with (3 fixed at j3 c . The order 
parameter should scale as (ipip) = dm 1 ' 5 . We accumulated data at quark masses ranging 
from 0.004 to 0.030, as shown in Figure ^| on a 12 3 x 6 lattice with the strength of the four 
fermi coupling set to 7 = 10. (We repeated this for 7 = 20.) Since the quark mass is different 
from zero in these simulations, the real order parameters (a), (ipip) can be measured rather 
than their absolute values. Our best fit gives a = 1.66(2) and critical index 5 = 3.89(3). 
The fact that this fit has only a 4% confidence is due to the point at m — 0.01 which lies 
2.5 standard deviations above the curve and is a problem for any smooth fit to the "data". 
Its removal increases the confidence to 33% with negligible change in the parameters of the 
fit. This lies between 3 and 5, the 5 values for the 2 leading symmetry breaking operators 
at the tricritical point. Since (ipip) could be expected to be a linear combination of these 
operators (and a third operator with 5 = 1 which is normally ignored), the more general 
scaling relation at (3 = (3 C is given by 

(^) = «0+^0 3 + ^c0 5 (2) 

with the solution of 

5 - m(a + b(f) 2 + c(f) 3 ) = (3) 

A fit to this form over the whole mass range yields a = 1.05(2), b = 2.5(4) and c = —3.0(1.4). 
Although the confidence level of this fit is only 2%, removing the m = 0.01 point improves 
this to an acceptible 34%. 
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FIG. 2. The order parameter plotted against m at criticality on a 12 3 x 6 lattice. The solid 
line is the fit (4>tp) = am 1 ^ . 

The hadonic screening lengths determine the manner in which chiral symmetry is 
restored as we pass through the transition from hadronic matter to a quark-gluon plasma. 
Here we concentrate on the 7r and cr(/ ) propagators which can be calculated from the tt and 
a auxiliary fields. For these we stored the ir and a fields averaged over each z-slice of the 
lattice, every 2 time units. Because the orientation of the condensate changes continuously 
throughout the runs, we relabeled the component of these fields in the direction of the 
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condensate as a and the orthogonal component as it for each configuration. These fields 
were then correlated to calculate the average a and tt screening propagators. The sum 
of these propagators is easily seen to give the joint a/n propagator in the plasma phase. 
The screening mass obtained from this unsubtracted propagator is plotted in Figure |3] as a 
function of (3. It flattens out in the hadronic matter phase when the fields develop a vacuum 
expectation value heralding a truly massless pion. While good fits of the pion propagator 
to a massless boson propagator are obtained in this region, fits in which the mass is allowed 
to float typically give masses ~ 0.2. To get clearer evidence for a massless pion would 
require a lattice of larger spatial extent (and/or better statistics). Until then we appeal 
to Goldstone's theorem for the knowledge that the pion is exactly massless in the chirally 
broken phase ||. The fit for the screening mass in the plasma phase is also shown. The 
fit used the data points at (3 = 5.4225 through (3 = 5.45 and used the form a{[3 — (3 C ) V ' . It 
predicted a correlation length exponent of v = 0.59(7), an amplitude a = 2.66(71) and a 
critical coupling j3 c = 5.4196(7). As is clear from the figure, the data here suffers from errors, 
both statistical and systematic, that are larger than in our other analyses. The prediction 
v = 0.59(7) is compatible with the 0.5 of the tricritical point. 

In summary, our measurements of (3 mag at N t = 6, made possible by our new action 
which permits zero mass simulations, strongly favour tricritical scaling at the chiral phase 
transition over the 0(4)/0(2) transition expected from universality arguements. Combined 
with the first order transition we have observed at JVj = 4 it gives hope that 0(4)/0(2) 
scaling will be observed for N t > 8. This 'anomalous' behaviour for N t = 4, 6 offers a 
potential explanation as to why simulations with the standard staggered fermion action 
produce unexpected results 0,[3],[1|] . 

The extra symmetry breaking operator at the tricritical point provides a possible 
explanation of the scaling of the chiral condensate with mass at (3 C and the measured value 
of 5. v is consistent with its tricritical value. 

Simulations on N t = 8 lattices should be performed to determine whether the contin- 
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FIG. 3. The pion screening mass plotted against (5 on a 18 3 x 6 lattice. The solid line is the fit 
discussed in the text. 

uum limit shows universal critical behaviour. Lattices with larger spatial extent are needed 
to obtain better predictions of v. 
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